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Abstract 

We introduce a model problem to address heat entrapment effects or the local accumulation of 
thermal energy within liquid acquisition devices. We show that the parametric space consists of six 
parameters, namely the Rayleigh and Prandtl numbers, the aspect ratio, and heat flux ratios for the 
bottom, side, and top boundaries of the enclosure. For the range of Ra considered 1 to 1 0 9 , beyond Ra on 
the order of 10 5 , convective instability is the dominant mode of convection in comparison to natural 
convection. The flow field transitions to asymmetric modes at Ra on the order of 10 7 . Direct numerical 
simulation of a large geometric length scale prototype for Ra on the order of 10 9 shows that the flow field 
evolves from small wavelength instability which gives rise to nonlinear growth of thermals, propagation 
of the instability occurs via growth of secondary and tertiary modes, and a travelling wave-type motion of 
convective modes occur prior to asymmetry. The effect of a large aspect ratio is to increase the number of 
modes in the vertical direction. Due to the slow diffusion of heat in the prototype, asymptotic states are 
not readily attained, we show that dynamical similarity can be used for a model which allows the 
attainment of asymptotic states and that transition to a chaotic state occurs for Ra on the order of 1 0 9 via a 
broadband power spectrum. These dynamical events show that for the baseline condition in which heat is 
absorbed from background laboratory environment, higher heat flux is absorbed at the top and bottom 
boundaries of the enclosure than a nominal value of 34.9 erg/cm 2 -s. 


Nomenclature 


Ra Rayleigh number 

Ar aspect ratio 

Pr Prandtl number 

f sidewall heat flux ratio 

fb bottom heat flux ratio 

ft top heat flux ratio 

q ” imposed heat flux 

H cavity height 

L cavity width 

g gravitational acceleration 

k thermal conductivity 

a thermal diffusivity 

p thermal expansion coefficient 

p density 

r boundary of cavity 
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Introduction 


Capillary liquid acquisition devices are used inside cryogenic storage tanks in order to supply vapor- 
free propellant to the engine. The liquid acquisition devise (LAD) interfaces directly with the feed system, 
for the supply of propellant, which can serve as a significant source of heat leak. In order to avoid 
excessive boil-off and ensure sufficient propellant subcooling to prevent cavitation and vapor formation 
during propellant outflow to the engine, the control of heat leak is necessary. Hence the localized 
accumulation of thermal energy, defined as heat entrapment, within the LAD flow channels becomes an 
issue of prime importance. We consider a model problem that addresses the effect of heat source at the 
boundary of an enclosure which contains typical cryogenic liquids (nitrogen, hydrogen) to lead insight 
into heat entrapment effects. This model problem has relevance towards the experimental works of 
Bolshinskiy et al. (Ref. 1), who considers heat entrapment effects modeled as controlled bottom heating 
of a cylindrical enclosure with and without partitions totally filled with a fluid. The partitions, which 
consist of typical screen meshes used in LADs as well as a solid divider, serve to emulate typical 
experimental conditions encountered during operations. Their results (Ref. 1) show that for the 
experimental conditions with partitions, the capillary screen meshes were impervious to natural 
convection. This is verified experimentally (Ref. 1) from comparison between data for meshes versus a 
solid boundary which show that the temperature rise are not significantly different. Hence a solid 
boundary can be used as an approximation of the boundary condition to address heat entrapment effects; 
this implies that the impermeable zero mass flux condition of the solid boundary models the boundary 
condition for the screen meshes. In our model we employ a solid boundary condition to address 
convection effects encountered in LAD applications. 

The basic physics of heat entrapment effects or the localized accumulation of thermal energy is 
addressed by broadly considering the effect of heat flux absorption on the boundary of an enclosure, due 
to controlled heat flux on the bottom and unknown heat leakage at the top and side boundaries, of an 
enclosed cryogenic liquid initially at saturation condition which totally fills the volume of the enclosure. 
Of basic interest is the rise in temperature due to various external heat flux conditions. The quantification 
of the temperature field can be used to assess optimum operating conditions for propellant flow to the 
engine. Related to the broad problem of the effect of heat absorption on the temperature rise of a 
cryogenic fluid are also issues on effect of thermal stratification on temperature rise due to side heating of 
partially filled cavities with a free surface as addressed by Tellep and Harper (Ref. 2), Vliet et al. (Ref. 3), 
and Vliet (Ref. 4) with applications to cryogenic storage tanks. Approximate boundary layer analysis 
using an integral approach coupled with experiments (Refs. 2 and 3) show that the effect of dominant side 
heating is to transport warmer fluid along the wall to the interface region which becomes accumulated at 
the free interface thus forming a stratified layer. Vliet (Ref. 4) addressed the effect of bottom heating 
using an approximate integral approach and show agreement with experimental data; the results show that 
as bottom heating increases there is a faster approach toward asymptotic conditions. Evans et al. (Ref. 5) 
have addressed the effect of non uniform heating on thermal stratification of water inside a rectangular 
enclosure and shown that thermal stratification can be reduced by distributing the heat flux along the side 
wall such that the maximum heat flux is concentrated near the bottom and decreases toward the top. The 
primary mechanism for the reduction of thermal stratification was postulated to be boundary separation at 
the wall. The bottom of the enclosure (Ref. 5) was insulated while the free interface was coated with a 
surfactant film of stearic acid to prevent vaporization; though the condition at the free surface is not 
specified, the effect of the surfactant film may be approximated as an insulating boundary condition based 
on the asymptotic profile of the rise in temperature of the enclosed fluid. 

The experimental works of Evans & Stefany (Ref. 6) addressed transient heat transfer effects inside 
totally filled cylindrical enclosures. They show that in response to a step-change in wall temperature, the 
mixed-mean temperature of the fluid exhibits an exponential rise in temperature which approaches the 
wall temperature asymptotically. In addition, the heat transfer coefficient oscillates in time as steady-state 
is approached. This finding was confirmed by the numerical and experimental works of Barakat & Clark 
(Ref. 7) for a partially filled cylindrical enclosure. Numerical works have been put forth by Polezhaev 
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(Refs. 8 and 9) who considers natural convection inside partially filled cavities for rectangular (Ref. 8) as 
well as cylindrical (Ref. 9) cavities for various boundary conditions. 

In contrast to the above works, our problem is somewhat unique in that there is no free surface and 
also that the fluid occupy the entire volume of the enclosure. The experimental works of Evans & Stefany 
(Ref. 6) share the similarity of total volume filling, however the boundary conditions are different, a step- 
change in temperature is prescribed. In our model we consider the Neumann prescribed heat flux 
conditions along the bottom, top and side of the enclosure. Thus depending on the dominant heat flux 
condition, there is coupling between natural convection driven by side heating and convective instability 
driven by bottom heating. The top heating serves as a stabilizing mechanism, due to imposed stable 
stratification, in contrast to bottom heating in which there is unstable stratification due to an increase of 
density with height. The dominant mechanism, natural convection or convective instability, depends on 
the external heat flux conditions during LAD operation. These two mechanism differ from the standpoint 
that in natural convection problems there is no stability threshold, convection occurs for any temperature 
variation across the boundary of the cavity whereas convective instability has a stable solution. The 
characteristics of natural convection inside vertical cavity due to differentially imposed temperature, at 
the vertical boundary, have been addressed by Elder (Refs. 10 and 1 1) for laminar (Ref. 10) and turbulent 
(Ref. 11) free convection. Transition from unicellular to secondary and tertiary flows have been shown to 
occur during laminar free convection, whereas the turbulent regime is characterized via travelling wave- 
like motion along the wall. The visual observation by Seki et al. (Ref. 12) has confirmed the findings of 
Elder, and also demonstrated the effect of Prandtl number. The mechanism of convective instability has 
been demonstrated for differentially heated horizontal cavity, with a forced flow component, whose 
bottom temperature is greater than the top by Akiyama et al. (Ref. 13). They show that for the limiting 
case of forced convection between horizontal plates, the limit of zero Reynolds number (no forced flow) 
yields agreement with the critical Rayleigh number criterion of 1708 corresponding to the classical 
Rayleigh-Benard problem (Ref. 14). Secondary flow patterns in the form of longitudinal vortex rolls were 
observed to occur near the critical Rayleigh number. The two basic configurations (Refs. 10 to 12) of a 
narrow vertical and horizontal cavity (Ref. 13) serve as basis by which the interaction between natural 
convection and convective instability mechanism can be understood. 

We introduce a model problem which considers the basic physics of heat entrapment inside an 
enclosure subjected to prescribed heat flux conditions to shed light on the convective state of the flow field 
in relation to the temperature field. This model problem is intended to lead insight into the temperature 
measurements (Ref. 1) made for liquid nitrogen inside an enclosure which absorbs heat from ambient 
laboratory condition, and also to assess the effect of imposed heat flux condition at the bottom boundary. 
Our main focus is to investigate the effect of heat flow from the ambient laboratory environment to the 
enclosure in order to establish a baseline to address the effect of imposed heat flux on the bottom boundary. 
The experiments were conducted using double encapsulation along the sidewall, thus the heat flux absorbed 
along the sidewall is minimal; we assume a value of heat flow on the order of 1 \lW. However, other 
unknowns are the heat absorbed on the top and bottom boundary for the experiments carried-out under 
ambient laboratory conditions. We investigate a range of heat flux condition to lead insight into its effects 
on the temperature rise of the fluid. This approach is necessary, since the baseline condition is not known 
and it is needed in order to assess properly the effect of bottom heating. 

We show that the effect of heat absorption is to destabilize the equilibrium condition of a fluid inside 
an enclosure and give rise to flow. The dynamical motion of the flow field obeys the evolutionary 
equations of mass, momentum and energy. For an incompressible Boussinesq fluid we numerically solve 
the coupled set of governing equations consisting of the Navier-Stokes coupled with the scalar 
temperature field using finite-difference methods. We consider the range of Rayleigh numbers (Ra) of 
order 1 to 1 0 9 for heating conditions comprising of prescribed heat flux at the sidewall, top, and bottom of 
the cavity. For the sidewall heat source we show that there is no stability threshold, incipient convection 
occurs. The flow field is dominated by two counter-rotating cells for low Rayleigh numbers and 
secondary flows can occur as Ra increases. However, for the bottom heating condition there is a stability 
threshold for which the onset of convection occurs which is on the order of 10 5 . As the Rayleigh number 
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increases, the horizontal wavelength which characterizes the onset of convection decreases. There is a 
transition to asymmetric flow for Ra on the order of 10 7 . Direct numerical simulation of the full scale 
prototype shows that for Ra on the order of 1 0 9 the onset of convective instability proceeds from short 
wavelength to nonlinear growth of thermals and culminates with a travelling wave-type mode prior to 
transition to asymmetry. Due to the large geometric length scale of the prototype, attainment of 
asymptotic states can be prohibitive. We employ dynamical similarity using a model with smaller 
geometric dimensions to show that similar dynamical events occur which facilitate the attainment of 
asymptotic states. As the Rayleigh number increases from 10 6 to 10 9 the flow field undergoes a transition 
from periodic to chaotic oscillations. The parametric study helps to identify the basic mechanism that 
drives the flow and lays the foundation to identify the baseline condition for heat entrapment effects. 

The paper is organized as follows: we postulate the basic physics of heat entrapment inside an 
enclosure. We derive the parametric space of the problem through scaling analysis. We solve the coupled 
set of nonlinear partial differential equations using finite difference methods coupled with a Flux- 
Corrected Transport method to resolve the small scale structure of the temperature field. We show that 
direct numerical simulation coupled with dynamical similarity can be used advantageously to model the 
large prototype enclosure used in the experiments. We show the basic convective state of the flow field 
for a relevant parametric space and that transition to asymmetric flow is relevant to this problem. 

Formulation 

The baseline condition for heat entrapment is shown fundamentally in Figure 1 A(a) and consists of a 
cryogenic liquid at saturation condition represented by p A which can absorb various amounts of heat at its 
boundaries comprising of surface heat fluxes q s , q t , qb denoting respectively the sidewall, top, and bottom 
boundaries. The enclosure in Figure 1 A can be approximated geometrically as the mid-plane of a cylinder 
with equivalent height (H) and diameter (L). This basic configuration represents the baseline condition 
investigated experimentally (Ref. 1) for liquid nitrogen inside a cylindrical Dewar under ambient laboratory 
condition; for which the heat flux absorb at the sidewall is controlled by external means such as double 
encapsulation of the wall region to minimize heat flow. Regardless of precautionary measures a finite 
amount of heat leakage occurred as indicated by the rise of temperature (Ref. 1) of liquid p A . This finding 
indicates that under laboratory environment, equilibrium condition is not attained for a given cryogenic 
liquid p A . One of the main focus of our model is to investigate the effect of relative magnitude of the heat 
fluxes (q s , q t , qb ) on the dynamical state of the liquid p A at saturation temperature T A when the system is 
driven off-equilibrium. This approach is necessary as it eliminates one of the unknown of the problem, if 
one of the heat flux is prescribed such as qb in the experiments (Ref. 1). Thus for a given qb , optimal values 
for q s and q t can be estimated based on ambient condition and experimental arrangement. Given the double 
encapsulation used on the sidewall of the experiments (Ref. 1), it is approximated that q s « q t . Hence for 
the baseline condition the dominant heat fluxes are qb and q t . 

More generally the effect of heat entrapment is investigated experimentally (Ref. 1) by placing a 
divider at mid-height, consisting of stainless-steel screen meshes or an aluminum solid divider, of the 
cavity as shown in Figure lA(b) in which p B = p A . The effect of heat entrapment is quantified 
experimentally by measuring the temperature above and below the divider which is used to compare 
against the baseline experiment illustrated in Figure lA(a). The findings (Ref. 1) indicate that 
temperatures below (bottom)/above (top) the divider position were higher/lower with the divider installed 
indicating heat entrapment/reduced stratification. The effects of the divider may be quantified by 
approximating the local heat flux q; into or out of the top region. Hence the problem may be simplified 
into two sub-problems as in Figure IB in which heat flows into or out of the top/bottom region shown in 
(a) and (b). This simplification effectively reduces the aspect ratio by 1/2 and the problem may be 
similarly addressed as the baseline problem in Figure 1 A(a) in which qb and q t can take the respective 
values of q; . This approach effectively reduces the parametric range of the problem which can be used to 
investigate the more general set of conditions as represented in Figure lA(b). 
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Figure 1 A. — Physical description of heat entrapment concept showing heat flux absorption (qb", q s ", qt") into 
cryogenic fluid (p A ) (a) without a solid divider and (b) with a solid divider which admits heat flux q" (-) into 
or out qj" (+) of the top region (p B ). 



Figure IB. — Subproblem of enclosure with solid divider showing (a) top region and (b) bottom region. 

In general, as shown in Figure 1 A(a), a condition for equilibrium to exist is for q s = q t = qb 
condition is not likely met for a low temperature cryogenic liquid inside an enclosure. Thus off- 
equilibrium condition need to satisfy the dynamical motion of the fluid govern by the continuity, 
momentum and energy equations as: 

du + dv _ q 
dx dy 


= 0; this 


( 1 ) 


du du du dp 

h u F v — = 1- v 

dt dx dy dx 


d 2 u 

_ 1 _ 

dx 2 

dy 2 ) 


dv dv dv 

h u 1- v = 

dt dx dy 
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( 4 ) 


df dT 8T d 2 T d 2 f 

— \-U— \-V~— = a 7T H 7T 

dt dx dy dx 2 dy 2 

In the above equations the density field is independent of pressure due to incompressibility of the liquid, 
however it is a function of temperature and can be linearized and given as an equation of state, 

p = p(r 4 )[i-(3Af] (5) 


We consider a Boussinesq fluid in which variation of density is significant only in the body force term as 
shown in Equation (3) where the coefficient of thermal expansivity is given as (3 = — . 

Superscript ~ denotes dimensional temperature and AT = T — T^. The thermophysical properties of the 
liquid, its density p, kinematic viscosity v, and thermal diffusivity a are taken to be at saturation condition 
for the specific cryogen, for liquid nitrogen T 4 = 77 °K. The gravitational acceleration g y is denoted by 
ng 0 in which n is a ratio by which the standard gravitational acceleration g G can be reduced via low 
gravitational environments such as the International Space Station (ISS), Moon, and Mars surface. 

The description of the initial and boundary conditions are made in reference to Figure 1 A(a). The 
initial condition is taken as uniform temperature throughout the liquid, 

t = 0,T(x,y,0) = T A ( 6 ) 


The momentum Equations (2) and (3) satisfy the no-slip condition along the boundary of the cavity 
denoted by T and given as, 


u{x , y,t) = 0 , v{x , y,t) = 0 on V 


(7) 


Whereas the energy Equation (4) satisfy the Neumann condition along the boundary F prescribed as, 


X = 0 ,q =qj ,0<y< H,9f/ dx)x=0 = - q"/k ( 8 ) 

x = L,q" =q" ,Q<y<H dT/ d3 ) x=L = -q" / k (9) 

y = 0, q = g&", o < x < L, d T/Q y ) y= o = -q b " /k (10) 

y = H, q = q t ", 0 < x < L, d T/^ y ) y=H = -q t " / k (11) 


Scaling Analysis and Computational Method 


The governing set of Equations (1) to (4) may be simplified by transforming the momentum 
Equations (2) and (3) into a vorticity equation using the following definition, 


u = 


d^jj 
8y ’ 




d'lfj 
dx ’ 


£ 


dv du 
dx dy 


( 12 ) 
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The resulting set of equations is scaled using the characteristics length, time, and velocity as 
(L,L / a, a / L ) , respectively. The dimensionless temperature is scaled as 


T — T. 

rp A_ 

q L / k 


The scaled dimensionless (denoted by *) set of equations become, 


aV aV 

dx* 2 dy 2 


( 13 ) 


(14) 


_ * 
dt 
_ * 

dt 


+ 


di * di 


dx 


dy 


= Pr 


'ay 

, 9y 

dx* 2 

dy* 2 


dT 

+ Ra Pr — *- 
dx 


dT 


dt 


* dT * dT 

U + V 


dx 


dy 


t q2 t q2 t 
+ 


[dx * 2 dy* 2 


The boundary conditions become, 

x =0,0 <y < Ar, &T/ ) * = — f 

/dx x =0 s 

x* =1, 0<y* <Ar = ~U 

y*=0,0<x*<l, d y dyt)y , =0 = -f b 

y* = Ar , 0 < x* <1,^^*)^ = -ft 

The set of parameters resulting from scaling is defined as, 

kva a L q q q 


(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 

(21) 


which represents respectively the Rayleigh and Prandtl numbers, aspect ratio, and dimensionless heat 
fluxes. Thus the parametric set, 


A = A(Ra,Ar,Pv,f s ,f b ,f t ) (22) 

consists of six parameters that can be varied independently and represents a formidable parametric set. 
Simplification of the parametric range is obtained by addressing the relevant limiting conditions for a 
given aspect ratio, fixed cryogenic liquid, and typical heat flux scale q . Since the predominant control 
variable is the bottom heat flux qb which implies fi,= 1, then f,f can take a range of values from zero to 
one. The conditions f = 0 and f= 1 signify insulated and heat absorption top boundary conditions, 
respectively. The side conditions f = 0 and f = 1 indicate a perfectly insulated sidewall and dominant 
sidewall heating. Since the sidewall has double encapsulation in the experiments (Ref. 1), we consider 
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values of f on the order of 1 0 6 which is approximately 1 p W of heat leakage. Although this heat leakage 
is insignificant in comparison to fab and fa it serves as a perturbation parameter to excite flow instability. 

The governing set of Equations (14) to (16) with boundary conditions (Eqs. (7) and (17) to (20)) are 
solved using finite-difference methods with a direct solver approach. The solution methodology consists 
of solving the stream function Poisson’s Equation (14) using matrix inversion which is followed by the 
solution of the energy Equation (16) and vorticity Equation (15). The vorticity field is used subsequently 
for each time step. Since we are interested in time evolution of the flow field dynamics, a time accurate 
methodology is employed which uses third order Adams Bashforth for time discretization. Since the 
bottom heat flux is dominant, the convective instability mechanism becomes important in comparison to 
natural convection, which requires high resolution of the temperature field. We employed the Flux 
Corrected Transport (FCT) method (Ref. 15) to resolve sharp interior gradients; in addition this method 
eliminates false numerical oscillation by ensuring positivity of the temperature scalar component 
throughout the domain of computation. 


Numerical Results 

Stabilization of Flow Field 

The dynamical state of fluid motion inside the enclosure shown in Figure lA(a) is first investigated 
for ideal cases in order to show the fundamental mechanism of convective motion. We first consider how 
stable states can be obtained and the cause of departures from stability. This is illustrated using liquid 
hydrogen as an example with Prandtl number Pr = 1.03. The destabilization of the flow field for bottom 
and side heat flux is summarized in Figure 2. The heating prescriptions that lead to convective stability in 
which there is no fluid motion even though conduction dominates is illustrated for Ra = 1.26xl0 3 in 
Figure 2(a). These trends show that stability of the flow field can be obtained, below the critical Rayleigh 
number, provided that there is a perfectly insulated sidewall fa = 0 for cases in which: (1) there is bottom 
heating only (fa= 1,/ = 0), (2) top heating only fa = 0 ,fi = 1), and (3) both top and bottom absorbing heat 
(f h = 1 ,fi= 1). The case when fa= \,fi= \) represents actual experimental condition, since neither top 
fa = 1 nor bottom heating fa = 1 can be isolated. However these cases represent idealization of an 
experimental system (Ref. 1) since there is finite heat leakage along the sidewall, i.e.,/> 0, regardless of 
precautionary measures. Note that absolute stability occurs for case 2 regardless of the Rayleigh number 
magnitude, however, at a critical Rayleigh number cases 1 and 3 become unstable marked by the onset of 
convection. 


Perturbation of Stability 

The perturbation of the stability for bottom fa= 1,/ = 0) and top heating (fb = 0 ,f = 1) is illustrated 
in Figure 2(b) for Ra = 1.26xl0 6 in order to isolate simultaneous effects of bottom and top heating, and to 
show the interaction of natural convection and convective instability. For a conservative heat leakage on 
the order of 1 p W projected on the sidewall surface of a cylinder fa = 2x1 0 -6 (qs = 0.001 erg/cm 2 -s), the 
stable base state fa= 1,/ = 0) becomes unstable due to unstable stratification caused by heavy fluid 
overlaying lighter fluid due to decrease of local density near the bottom of the enclosure caused by 
bottom heating (fa = !)• This case corresponds to qb =500 erg/cm 2 -s which can be obtained from heat 
flow from ambient laboratory environment on the order of 1 0 mW. In this case convective instability 
driven by bottom heating is dominant over natural convection driven by side heating. Since we are 
operating near the threshold of instability the growth of four modes occur. In contrast to instability of the 
bottom stratification fa= l,fa= 0 ,fi = 0), with a side heat leakage of fa = 2x1 CT 6 , we show that the top 
stratified layer fa = 0,f = 1,/ = 0) remains stable for the same heat leakage; however horizontal density 
gradient gives rise to a natural convection base flow with a symmetric double cell. These two cases 
illustrate that even though the top stratified layer is stable for all Rayleigh numbers when the sidewall is 
insulated/* = 0, the effect of finite side heat leakage is to drive a natural convection base flow. 
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Figure 2. — Destabilization of flow field due bottom and side heat flux (t = 360 s); qb", qt" are nominally 500 erg/cm 2 -s, 
q s " = 0.001 erg/cm 2 -s for side heat flux perturbation. 
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Effect of Side Heating 


Since side heat leakage is an effective source to excite instability and drive a natural convective base 
flow, we isolate the effect of side heating by considering the case when there is no heat leakage from the 
bottom (fb — 0) and top (f, = 0) contrasting Ra = 2.5 and Ra = 1 .26x 10 6 in Figure 2(c). These two cases 
contrast heat flow on the order of 1 pfT(qs” = 0.001 erg/cm 2 -s) and 10 mJT(qs” = 500 erg/cm 2 -s) 
corresponding to Ra = 2.5 and Ra = 1.26xl0 6 , respectively. For Ra = 2.5 natural convection becomes 
dominant, it generates a very weak flow field which does not disturb the conduction layer near the vertical 
walls. However, when there is significant side heat leakage, Ra = 1.26xl0 6 , the cellular structure becomes 
unstable which results in growth of impending secondary modes due to shear instability. There is a 
transition from conduction dominated regime to a boundary layer regime indicated by the temperature field. 
The flow field transports warm fluid near the walls to the interface region which forms a stratified layer. The 
scenario of sidewall dominated heating is similar to events that occur in applications to address the effect of 
thermal stratification for partially filled cryogenic storage tanks (Refs. 8 and 9). The case for Ra = 1.26xl0 6 
(fb = 0,f = 0) is demonstrated for illustrative purposes; the precautionary measures taken to insulate the 
sidewall in the experiments (Ref. 1) would indicate a much smaller side heat leakage toward the lower limit 
of Ra = 2.5. 


Perturbation of Double Layer 

A typical heating condition with finite sidewall heating is illustrated in Figure 2(d) for the case when 
there is simultaneous heat absorption from the top and bottom boundaries which represents perturbation 
of the double layer for Ra = 1.26xl0 6 based on case 3. The destabilization of the conduction dominated 
regime (fi,= 1 ,f= l,f s = 0), near the critical Rayleigh number, is contrasted for/ 9 = 2x1 CT 6 and f s = 1 . 
When f = 2XKT 6 convective instability dominates, there occurs growth of four modes at the bottom 
region while the growth of the stratified layer at the top region for early time t = 360 s remains stable. 
However, if there is significant sidewall heat leakage (fb= l,f - 1 ,f= 1) natural convection dominates 
and the characteristic of convection switches from an instability mechanism to a boundary layer 
dominated flow shown by the convective cells along the wall and secondary cells near the bottom. In this 
case the intensity of convection modifies the top stratified layer and there is a potential for local 
superheating near the top comers. In summary, the mechanism by which the double layer (fb = 1 ,f= 1, 
f s = 0) becomes unstable is illustrated by perturbation of the stability with f = 2x1 CT 6 , and the effect of 
dominant side heating on the stability is illustrated for^ = 1 . This indicates that there can be a wide range 
of dynamical events that occurs for 2x10 ~ 6 <f s < 1, even though we fix/* = 2x 10~ 6 for the parametric 
study, it is a key parameter to characterize heat entrapment effects. The results also show that when the 
double layer is perturbed convective mixing becomes important. These limiting conditions serve to 
illustrate the interaction of convective instability and natural convection that can occur for the baseline 
condition shown in Figure lA(a). 


Mechanism of Instability 

For the given experimental conditions, it is likely that the convective characteristic inside the 
enclosure is dominated by convective instability (f = 2XKT 6 ). In Figure 3 A we illustrate the mechanism 
of the instability as the Rayleigh number increases. The results show that from the threshold of instability 
Ra = 1.26xl0 5 there is a transition from long to short wavelength as Ra increases to 1.26xl0 9 ; there is 
saturation in mode numbers between Ra = 1.26xl0 7 and Ra = 1.26xl0 9 . The growth of modes segregates 
near the walls as Ra increases and the wavelength between the cells decreases. This trend indicates that 
high Rayleigh number flows require resolution of small scales for computational accuracy. 
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Ra = 1.26x105 
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Ra= 1.26x109 
T max = 0.05, t = 12 s 
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Figure 3A. — Transition to short wavelength instability as Ra increases 
(f b = 1 , f t = 1 , f s = 2x10-6) q" = 500 erg/cm2 -s, Ar = 1 , Pr = 1 .03. 
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Ra 


Figure 3B. — Horizontal wavelength of flow field as a function of Ra (fb = 1 , ft = 1 , 
f s = 2.0*10-6) Ar = 1, Pr= 1.03. 

The characteristic wavelength (X = X * / L) as a function of the Rayleigh number in Figure 3B 
shows that X is inversely proportional to Ra, thus decreases sharply as Ra increases. For the highest 
Rayleigh number considered (1.26xl0 9 ), A = 0.075 which is on the order of 4 mm. This finding indicates 
that the length scale which needs to be resolved gets smaller and smaller for Ra > 10 9 and also illustrates 
the challenge posed for computations of high Rayleigh numbers. 

Dynamical Events 

We now turn to the dynamical events as the flow field approaches asymptotic states. This is 
illustrated for Ra = 1.26xl0 6 in Figure 4, which shows growth of a single thermal, and Ra = 1.26xl0 7 in 
Figure 5 A which shows growth of double thermals and the breakdown of symmetry in Figure 5B. In 
contrast to natural convection dominated flows in which the maximum temperature occurs near the 
sidewalls due to side heating (Fig. 2(c)), in this convective instability dominated situation (Fig. 4), the 
maximum temperature occurs near the mid-center of the cavity due to bottom heating since the sidewalls 
are nearly insulated (f = 2x1 0 -5 ). The mechanism of the instability is due to the release of potential 
energy from unstable stratification caused by heavy fluid overlaying lighter fluid and its conversion to 
kinetic energy to drive fluid motion. In this case, the warmer fluid transported upwards and cooler fluid 
moving downwards result in the release of potential energy; the inner cells grow faster than the outer cells 
since the inner region has the maximum temperature. The inner cells grow as time increases, thus the 
instability self propagates until the top wall is reached. The initial bottom four cells (t = 336 s) self- 
organizes into two rows of secondary cells ( t = 408 s) due to the nonlinear growth of the thermal. Since 
the top of the cavity is insulated (f t = 0), the intensity of convection increases as the top of the cavity is 
approached. Buoyancy effects due to growth of the thermal propel the two inner cells toward the top of 
the cavity (t = 360 s) while the outer cells near the walls segregate toward the center of the cavity 
(t = 384 s). The maximum temperature decreases with time due to local mixing induced by the increase in 
convective intensity for Ra = 1.26xl0 6 . 
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Figure 5A. — Evolution of thermals due to nonlinec 
Ra = 1.26x107, Ar= 1, Pr= 1.03. 
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Figure 5B— Transition to asymmetric flow, (fb = 1, ft = 1, fs = 2x1 Ch 6 ), Ra = 1.26x10 7 , Ar = 1, Pr = 1.03. 

We consider the dynamical evolution of the flow field when there are six modes in which heat 
absorption occurs at the top boundary as shown in Figures 5 A and B for Ra = 1.26xl0 7 . There is 
similarity of events in comparison to the top insulated case for Ra = 1.26x1 0 6 . Since the maximum 
temperature occurs toward the mid-section of the cavity, the inner cells grow faster than the outer cells, 
the growth of six modes leads to formation of two symmetric thermals with respect to the center-line. The 
instability produces buoyancy which lifts the two central cells thus forming tertiary modes (t= 130 s); the 
rows of cells from the center toward the wall form partitioned co-rotating cells with the middle cell 
(tertiary mode) counter-rotating. The tertiary modes propagate toward the top of the cavity (/ = 135 s) 
from a bifurcation of a local double homoclinic orbit which signals transient effects. There is subsequent 
merging of the smaller cells (tertiary modes) along the bottom, forming a single homoclinic orbit, into 
single cells (/ = 140 s). The growth of the stratified layer from the top limits the penetration of the modes 
as shown in Figure 5B (/ = 960 s); thus the system transitions to four modes. These four modes become 
asymmetric ( t = 1080 s) and the growth of the asymmetry leads to merging of the inner and outer cells 
(secondary modes). Thus two asymmetric modes are formed at t = 1200 s. This finding is significant as it 
implies that beyond the critical Rayleigh number the flow field becomes asymmetric. This implies that 
axi-symmetric models would be limited in their application to high Rayleigh numbers since this occurs 
for Ra = 1.26xl0 7 . 
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Local Measures 


The local measures which characterize the behavior of the fluid as a function of system parameters 
are summarized in Figures 6 A, B, and C. The effect of increasing the Rayleigh number is considered in 
Figures 6 A and B by showing the characteristic maximum velocity magnitude and temperature as a 
function of Ra contrasting the effect of insulated (f = 0) and top absorbing (f= 1) boundaries at t = 600 s. 
The results in Figure 6 A show that V max remains constant for both fi = 1 and f — 0 for Ra < 10 5 indicating 
stability. Instability occurs for Ra > 10 5 for which the solution remains the same up to 10 7 . There is faster 
rise in velocity for the case f = 0 for 10 7 < Ra < 10 9 . This represents the range for 



Ra 


Figure 6A. — Characteristic velocity magnitude at t = 600 s, comparing heat flux absorption (fy = 1 ) 
and insulated (ft = 0) top boundary conditions for f b = 1 , f s = 2. Ox 1 0 -6 , Ar = 1 , Pr = 1 .03. 



10 ° 10 1 1 0 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 

Ra 

Figure 6B. — Characteristic maximum temperature at t = 600 s, comparing (ft = 1 and 
ft = 0) top boundary conditions for fb = 1, f s = 2.0x1 0~ 6 , Ar = 1, Pr = 1.03. 
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Figure 6C. — Transient response of cryogenic fluid comparing heat absorption (ft = 1) to 
insulating (f t = 0) top, f b = 1, f s = 2.0x10-6, Ra = 1.26x106, Ar = 1, Pr = 1 .03. 


which the stabilizing effect of the top boundary / = 1 becomes effective and lowers the convective velocity. 
The trend of the maximum temperature as a function of Ra in Figure 6B shows that T max is higher for/ = 0 
when Ra < 10 7 whereas for Ra > 10 7 T max becomes greater for the stabilizing top heat absorption / = 1 case. 
These trends indicate that local convective mixing is much more effective at lowering T max for the insulated 
top as opposed to the case when there is heat absorption on top for Ra > 10 7 . In terms of heat entrapment 
(localized increase in temperature), the trend from / = 1 to/ = 0 is representative of expectations, that is an 
increase in local maximum temperature for Ra < 10 7 and /< 1; since the effect of the divider is to limit heat 
flow, the region bounded by 0 </< 1 represents potential conditions where T max increases beyond the 
limiting condition / = 1. The effect of the divider approximates the trend for/< 1, thus / = 0 represents a 
limiting condition for heat entrapment. For Ra > 10 7 the increase in convective effects for/ = 0 induces 
local mixing which decreases heat entrapment effects, thus a reverse trend is predicted with an asymptotic 
approach to a lower temperature in comparison to/ = 1. If heat entrapment is based on the condition /< 1, 
indicating finite heat leakage at the divider boundary then heat entrapment is less effective for Ra < 10 7 
when / = 1, whereas for Ra > 10 7 it becomes more effective. The trend for Ra > 10 7 is in agreement with the 
trend in experimental measurements of temperature increase near the bottom divider (Ref. 1), since Ra > 10 7 
for the experimental system. In Figure 6C, we show the transient response of a cryogenic fluid along the 
central axis (0.5, y) of the cavity for/ = 1 and / = 0. The trend of an adverse temperature gradient at the 
bottom of the enclosure is not greatly affected by the top conditions. The effect of top heating (f= 1) is to 
create a stably stratified layer which increases with time. Whereas, when the top is insulated / = 0, there is 
no stratified layer; the temperature becomes nearly uniform near the top as the asymptotic state is 
approached. 


Direct Numerical Simulation/Dynamical Similarity 

We demonstrated the characteristic of the flow field and the response of the temperature field for a 
fixed Pr, Ar, and a range of Rayleigh number for a model. We now consider the prototype used in the 
experiments (Ref. 1) which has much greater geometric length scales, that is L = 19.1 cm, H = 61.9 cm. 
Due to the large length scale employed in the prototype it requires a much larger grid density to resolve 
the small length scales in the flow and temperature fields that occurs for large Rayleigh numbers Ra > 10 8 
as shown in Figure 3B. For example the horizontal length scale of the prototype can be four times greater 
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than that of the model; for a typical grid size used for the model of 200x200, the prototype would require 
a grid density four times as large of 800x800 to resolve the same length scale as the model. This can 
require more than doubling the computational resources; a typical case for a model which requires 
between 30 to 72 h using 4 Central Processing Units (CPU’s) now requires between 300 to 800 h using 
20 to 32 CPU’s for the prototype. Even though, direct numerical simulation of the prototype is the most 
exact, dynamical similarity can be used to model the prototype. This leads to considerable savings in 
computational resources. We explore the application of dynamical similarity by showing that a model 
with smaller geometric length scale, requires matching of the parametric set A p (given by Eq. (22)) for the 
prototype to that of the model A M ; we also show that dynamical events are similar with the exception of 
the time scale. The issue being that due to reduction of length scale of the model, the equivalent heat flux 
diffusing through the model takes much shorter time. The model allows the dynamical events to be 
accelerated, hence the savings in computational resources. 

In order to address direct numerical simulation of the prototype (for a geometric length scale 
L = 19. 1 cm, H = 61 .9 cm), we consider Ar = 1 as a base case, Ar = 1 .32 (H/2) for the half scale model in 
Figure IB, and Ar = 3.26 (H) for the full length of the prototype in Figure 1A. The external condition for 
the prototype is absorption of heat on the order of 1 mW for the top and bottom boundaries, due to 
ambient background laboratory condition, which is equivalent to its projection on the surface of a circle 
with heat flux of q = 34.9 erg/cm 2 -s for liquid nitrogen. The condition at the sidewall is heat absorption 
on the order of 1 p W or a heat flux q” = 0.001 erg/cm 2 -s. The characteristics of the dynamical evolution 
of the flow field is illustrated for Ar = 1 in Figures 7 A, B, and C which show three main events: (a) 
transition to short wavelength instability and growth of thermals, (b) vertical propagation of the 
instability, and (c) a travelling wave-type mode which eventually transitions to asymmetry. 

Even though the cavity size for the prototype in Figure 7 A is much larger than a typical model shown 
in Figure 3 A (Ra = 1.26xl0 9 ), for nearly equivalent Rayleigh numbers and slightly different Prandtl 
number, similar transition instability occurs in Figure 3A at t = 12 s as compared to Figure 7A at t = 288 s 
for Ra = l.OxlO 9 and Pr = 2.27. In this case there is a deceleration of the instability since the size of the 
cavity is so large for a given heat flux input. There is a similar growth of thermals to that shown in 
Figure 5A, however the higher grid density used (500x500) allows the resolution of small length scales 
such as the fine structures during the interaction between thermal pairs shown in Figure 7 A at t = 432 s. 
The propagation of the instability toward the top of the cavity occurs via a repetitive sequence of growth 
of secondary cells in Figure 7B. The mechanism of the instability occurs through the formation of 
homoclinic orbits ( t = 648 s and t = 720 s), which serve to increase the size of the cells and form 
secondary cells inside the main cell, which subsequently detaches and form secondary cells as shown 
from the transition between t = 720 s to t = 792 s. Once the flow field penetrates near the top of the cavity, 
there is a travelling wave-type mode event that occurs as shown in Figure 1C. The travelling wave-type 
mode occurs through the exchange of energy between tertiary and secondary modes as shown for 
t = 864 s, t = 1296 s, and t = 1368 s. As the travelling wave-type mode propagates up and down the cavity 
there is a breakdown and reconnection of the vertical modes; there is also boundary layer entrainment 
along the vertical walls. The flow field subsequently transitions to asymmetry. Note in comparison to 
Figures 5 A and B the boundary layer extends further up the cavity vertical walls. 

Effect of Aspect Ratio 

The effect of increasing the height of the cavity to actual prototype dimension model Ar = 1.62 (H/2) 
and Ar = 3.24 (H) for long and short time is shown in Figure 8 contrasting insulated/ = 0 and absorbing 
top f = 1 boundary conditions. For the case of a taller cavity Ar = 1.62, the growth of tertiary modes is 
formed between two co-rotating secondary cells. The tertiary mode has a homoclinic orbit which signals 
propagation of a travelling wave inside a fluid. This case is contrasted with Ar = 3.72, the prototype used 
in the experiments (Ref. 1). Such a large cavity requires high grid density to resolve small scales; we 
show direct numerical simulation of early time evolution of the instability using 1000x1000 grids for 
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short times. The results show the onset of instability which consists of six modes segregated toward the 
wall with small length scales. In this case, the characteristic time for the instability to occur increases 
(t = 288 s) due to the large increase in vertical length scale as compared to a model in Figure 3A which 
shows the instability at t = 12 s; note the instability develops faster with a lower aspect ratio prototype as 
in Figure 7A at t = 288 s. 


T(x,y,t) T(x,y,t) ||V||(x,y,t) 

Tmax = 0.02, V max = 1 .9x10 -6 cm/s, t = 72 s 



Tmax = 0.47, V max = 0.05 cm/s, t = 288 s 

1.0 T T T T 

0.8 - 
0.6 
0.4 


0.2 





Figure 7A. — Transition from short wavelength instability to nonlinear growth of thermals (fb = 1 , ft = 1 , 
f s = 2.9x1 0~ 5 ) for prototype (L = H = 19.1 cm); q" = 34.9 erg/cm 2 -s, Ra = 1.0x10 s , Ar = 1, Pr = 2.27 
(grid size 500x500). 
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Figure 7C. — Travelling wave-type motion of convective modes, Ra = I.OxlO 9 , Ar = 1, Pr = 2.27. 
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Figure 8. — Effect of aspect ratio on propagation of instability, (Ar = 1 .62, H = 30.95 cm) and onset of 
short wavelength instability (Ar = 3.24, H = 61.9 cm) for prototype (L = 19.1 cm), q" = 34.9 erg/cm 2 -s, 

Pr = 2.27, Ra = 1.0x109. 

Application of Dynamical Similarity 

We illustrate the application of dynamical similarity using a model in Figure 9 with smaller length 
scale for the same heat flux boundary condition as the prototype in Figure 7 A. In this case, 

A M (Ra,Ar,PT,f s J b J t ) = A p (Ra,Ar,PYj s J b ,f t ) (23) 

note the matching of the parametric space yields similar dynamical events as compared to Figures 7A, B, 
C except with the model the characteristic time for the instability to develop is much shorter t = 20 s as 
compared to / = 288 s for the prototype in Figure 7 A. 

Local Flow Dynamics at Fixed Points 

The local dynamics of the flow and temperature fields at fixed points (x p , y p ) as denoted in 
Figure lA(b) is shown in Figure 10. The decomposition of the time history signal into its frequency 
components for the flow field is obtained by estimating its power spectrum using, 

T h 2 

W) = Jr \v(i)e~ 2vift (24) 

T b 0 

in which V(t) is the time history of the corresponding component of the velocity field at the time interval 

of the data denoted by Tb. The smooth estimate Pv(t) is obtained from a convolution relationship 

(. P v (/) = W H (/) * Py (/)) using a Hanning spectral window (W H (f)). The dynamical characteristics of 
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Figure 10. — Dynamical characteristics of flow and temperature fields at various fixed points (x p ,y p ) for prototype, 
Ar =1,Pr = 2.27, Ra = 1.0x109 (grid size 500*500). 
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flow and temperature fields at various fixed points (x p , yp) for the prototype are shown in Figure 10. In 
relation to Figure 1 A(b), the series of points represent location near the left wall which proceeds from the 
bottom toward the top for a time interval Tb = 1500 s. As shown in Figure 10, the time history of the 
temperature and flow fields in the bottom region at (0.002, 0.002) shows that for a typical q” = 

34.9 erg/cm 2 -s in which (ft, = 1 ,f= 1 ,f = 2.9x1 0 -5 ), the temperature increases initially by conduction 
effects up to t = 300 s while the flow field remains quiescent; for t > 300 s the onset of instability occurs 
in which the flow and temperature fields oscillate aperiodically in the core region. As the distance 
increase from (0.15, 0.1) to (0.15, 0.35) at the bottom, it takes longer ( t = 400 s) for the temperature and 
flow fields to respond and there is also an increase in their aperiodicity. As the distance is further 
increased toward the top of the cavity, from (0.15, 0.65) to (0.002, 0.998), the flow and temperature 
oscillations transition from aperiodic to a nearly periodic (mono-frequency) oscillation. This trend is due 
to absorption of heat (f t = 1) at the top boundary which stabilizes the flow. The power spectrum shows 
that there is a dominant mono-frequency component of 0.004 Hz with slight harmonics at (0.15, 0.90); 
this corresponds to a long period of 250 s. Note that near the top (0.002, 0.998) conduction effects 
dominate and it takes on the order of t = 800 s for the onset of convective instability. The behavior of the 
wall temperature at the mid-section of the bottom (0.5, 0.0) and top (0.5, 1.0) shows that at the bottom 
convective instability ( t > 400 s) results in a decrease in temperature due to cooling by the flow field 
whereas at the top the temperature increases beyond the onset of oscillations ( t > 800 s). These trends are 
representative of the behavior of the flow and temperature fields for the prototype for Ra = l.OxlO 9 . 

We now consider the behavior of the flow and temperature fields for a model near the threshold of 
instability at Ra = 1.26xl0 6 as well as a case for Ra = l.OxlO 9 in order to compare with the prototype. The 
application of dynamical similarity (Am = A P ) using a model with smaller geometric length scale, showing 
time asymptotic behavior near the top of the cavity at (0.15, 0.90), is shown in Figure 11. Since the cavity 
size of the model is smaller, the equivalent amount of heat input to the bottom diffuses much faster as 
compared to the prototype, therefore the instability occurs much sooner. Near the threshold of instability 
at Ra = 1.26xl0 6 there is oscillation in both the temperature and flow fields; the power spectrum of the 
flow field shows oscillations with frequency of '/= 0.001 Hz and harmonics of the frequency, 2 f 3 f 4 f 
(4/= 0.004 Hz) which spans a range of period from 250 to 1000 s. Comparison of the model, showing a 
time asymptotic state, and prototype in Figure 10 at location (0.15, 0.9) shows that the model exhibits 
aperiodic behavior as time increases for Ra = l.OxlO 9 . The time asymptotic behavior of the velocity 
component shows intermittent periodic bursts which contain a broadband frequency content extending to 
0.4 Hz as indicated by the power spectrum P u ; this dynamical state indicates a transition to chaos. We 
compare a sub-interval of the time history to the prototype in Figure 10 at (0.15, 0.90) and the remarkable 
results show similar behavior for both the temperature and flow fields. The initial response time of the 
model occurs near t = 10 s whereas the prototype corresponds to t = 200 s, and the onset of convective 
instability occurs near t = 60 s whereas the prototype indicates t = 700 s. Similar to the prototype the 
power spectrum indicates periodic frequency response, however the frequency (0.04 Hz) is ten times 
larger with a corresponding shorter period of 25 s, and higher harmonics indicating nonlinear effects. 

Such dynamical similar results indicate that a model can be used to infer the dynamical characteristics of 
a prototype approximately. This finding is practical as it represents an alternate means to study the 
behavior of a prototype which requires significantly less computational resources; for example a typical 
prototype that requires 800 h to run may only require 40 h for a model. The major difference between the 
model and prototype is the characteristic time for the onset of instability, the model requires much shorter 
time, hence the savings in computational time. 
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Figure 1 1 . — Similarity in dynamical characteristics for model (Ar = 1 , L = H = 5 cm) at various Rayleigh numbers, 
(grid size 200x200). 


The effect of aspect ratio on the prototype for Ar = 1 .6 and Ar = 3.24 is contrasted for long and short 
times in Figure 12. For Ar = 1.6, we show a case for Ra = 1.0x1 0 9 corresponding to subcooling at the top 
boundary (f = — 0.1) which is relatively low. The time history of the temperature and flow fields near the 
bottom of the enclosure at (0.0025, 0.0017) is similar to that in Figure 10 shown at (0.002, 0.002). Thus low 
subcooling on the top boundary does not have a significant effect on the dynamical characteristics of the 
flow field. We show a full scale simulation of the prototype for Ar = 3.24 in which the grid density is 
increased from 400x600 (Ar = 1.6) to 1000x1000; this is demonstrated for the baseline condition in which 
the same amount of heat flux is absorbed at both the bottom and top boundary (fb = l,f= 1) for short times. 
The results show that it takes on the order of 300 s for the onset of instability; this corresponds to 
approximately 800 h of computational time. To compute the time asymptotic behavior requires significant 
resources to reach the asymptotic state. The wall temperature at the bottom (0.5, 0.0) shows a monotonic 
increase without fluctuations; this occurs because for the short time considered conduction effects are 
dominant. In comparison to the parametric range of the experiment (Ref. 1), the upper limit of imposed heat 
flow at the bottom boundary is 104 W which corresponds to approximately qb = 3.49xl0 6 erg/cm 2 -s. In 
order to find the baseline condition to estimate the heat flux absorbed from the top and bottom boundary we 
assumed a heat flow of approximately 1 mW, from background environment, which corresponds to a heat 
flux of qb = 3.49x10' erg/cm 2 -s. The scaling of the temperature magnitude using Equation (13), for 
comparison with experimental measurements (Ref. 1) for the baseline condition of no screen, indicates that 
the heat flux absorbed at the boundary is greater than the assumed value of qb = 3 .49x10* erg/cm 2 -s and lies 
in the range of 3.49x10* erg/cm 2 -s < qb < 3.49x1 0 6 erg/cm 2 -s. Using a combination of dynamical 
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Ar = 1.6, (0.0025, 0.0017) (400><600) 



Ar = 3.24, (0.0025, 0.0017) (1000x1000) 
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Figure 12. — Effect of aspect ratio on dynamical characteristics of prototype with (Ar = 1.6, L = 19.1 cm, 
H = 30.95 cm) and (Ar = 3.24, L = 19.1 cm, H = 61 .9 cm) for Pr = 2.27, Ra = 1.0x109. 


similarity and full scale simulation of the prototype, it remains to investigate the effect of higher heat flux 
absorption at the boundaries. However the parametric range considered, for low heat flux condition at the 
top and bottom boundaries, provides a foundation to investigate heat entrapment effects at higher Rayleigh 
numbers. 


Summary and Conclusions 

We address heat entrapment effects with application to liquid acquisition device (LAD) inside 
cryogenic storage tanks. We use dynamical similarity for the model along with direct numerical 
simulation of the prototype with much larger geometric length scale to address the global dynamics of the 
flow and temperature fields and local dynamics at fixed points. From scaling analysis, we show that the 
parametric space depends on six parameters namely (Ra, Pr, Ar ,f s ,fb,ft) the Rayleigh and Prandtl 
numbers, the aspect ratio, and the corresponding heat flux ratios for the side, bottom, and top boundaries. 
Due to the large parametric set involved, we show a selected parametric range based on typical 
experimental conditions applicable to LAD. The parametric set entails Ra between 1 and 1.26xl0 9 , Ar 
between 1 and 3.24, Pr = 1.03 and 2.27, (f s ,fb ) between 0 and 1, and f between -0.1 to 1. We address 
particularly the effect of applied heat flux at the boundaries to establish a baseline condition for a non- 
equilibrium cryogenic fluid inside an enclosure absorbing heat due to background laboratory 
environment. This serves to establish the conditions likely to occur at the top and side boundaries which 
can be used for simulation with known input heat flux to the bottom boundary. 
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The results show that for typical conditions, convection is driven via the interaction between natural 
convection driven by heat flux absorbed along the vertical walls and convective instability due to heat 
flux absorbed along the horizontal bottom boundary. Heat flux absorbed at the top boundary serves as a 
stabilizing mechanism. The predominant mechanism for convection is convective instability for which the 
longest wavelength occurs for Ra on the order of 10 5 which is near the threshold of instability. The 
maximum number of modes that occur is six for Ra of order 10 7 ; these modes segregate toward the wall 
as Ra increases to 1.26xl0 9 . The wavelength of the horizontal modes decreases as Ra increases which 
governs the length scale that needs to be resolved for computations of high Rayleigh numbers. For 
Rayleigh numbers near the threshold of instability the flow field remains symmetric, however for Ra on 
the order 10 7 there is a transition to asymmetric flow. 

Direct numerical simulation for large geometric length scale prototype for Ra = l.OxlO 9 shows that 
the global dynamics of the flow field evolve in three stages. These stages consist of a transition from short 
wavelength instability to nonlinear growth of thermals, propagation of the flow field to the top of the 
cavity via a repetitive sequence of destabilization which gives rise to secondary and tertiary modes, and a 
travelling wave-type motion of convective modes which occur prior to transition to asymmetry. The 
effect of increasing the aspect ratio is to increase the number of modes in the vertical direction. Due to the 
slow diffusion of heat for large aspect ratio enclosures, the onset of instability occurs much later in 
comparison to a model with smaller geometric dimensions and computation of asymptotic dynamics of 
the flow field becomes challenging for the prototype. We use dynamical similarity to show that similar 
events can be obtained with a model which facilitates the attainment of asymptotic states. The time 
history at specific points in the flow field shows that for Ra on the order of 10 6 the temperature and flow 
fields oscillate periodically, however for Ra on the order of 10 9 there is aperiodic behavior with 
intermittency in the flow field which is indicative of a transition to a chaotic state. Scaling of the 
temperature magnitude, for the parametric range of the model and prototype, indicates that there is a 
higher heat flux absorption on the boundaries than the nominal value of 34.9 erg/cm 2 -s or approximately 
1 mW of background ambient heat flow. However, the parametic study shed light on the basic flow field 
dynamics associated with heat entrapment and can be used in future studies to address effect of increased 
heat flux at the boundaries. One of the highlights of the parametric study is the finding that transition to 
asymmetry occurs for Rayleigh numbers relevant to heat entrapment, this implies that axi-symmetric 
models are to be used with caution, since they are applicable to a very narrow parametric range near the 
critical Rayleigh number. 
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